// -* Data Analysis for Belle Experiment -*-
//
// Description: B0 -> gamma gamma 
//
// Author : Yuuji Unno, Taeyoung Kim
// Created: 2009/04/07.
//
//
// Revision History
//
// Revision 0.0  2009/04/07      Y.Unno
//    First building the Program
//
// Revision 0.1  2009/04/10      T.Y.Kim
//    add Momentum, Energy
//
// Revision 0.2  2009/04/16      T.Y.Kim
//    add theta angle
//
// Revision 0.3  2009/05/01      T.Y.Kim
//    add phi angle
//
// Revision 0.4  2009/07/02      T.Y.Kim
//    add e9/e25
//
// Revision 0.5  2009/07/16      T.Y.Kim
//    add timing of TC & csi_bb
//
// Revision 0.51 2009/09/02      Y.Unno
//    add eventfunction, rec_b2gg
//
// Revision 2.0  2009/09/04      T.Y.Kim
//    add SCF
//
// Revision 3.0  2009/09/17      T.Y.Kim
//    add Direct pi0,eta veto
//
// Revision 3.1  2009/09/25      T.Y.Kim
//    add pi0 using mdst_pi0
//
// Revision 3.11 2009/10/10      T.Y.Kim
//    compare pi0 using mdst_gamma, mdst_pi0
//
// Revision 3.12 2009/10/17      T.Y.Kim
//    add # of candidate
//
// Revision 4.0  2009/11/06      T.Y.Kim
//    add helicity angle, asymmetry gamma energy
//
// Revision 4.2  2009/11/13      T.Y.Kim
//    del helicity angle, asymmetry gamma energy
//    add eta mass
//
// Revision 4.3  2010/01/20      T.Y.Kim
//    del pi0 using mdst_pi0
//    add mc_gg
//    change pi0,eta -> gg
//
// Revision 4.4  2010/01/31      T.Y.Kim
//    add pi0etaveto function
//    
// Revision 4.5  2010/02/14      T.Y.Kim  
//    add pi0etaveto, applying 2 gamma
//
// Revision 4.6  2010/03/06      T.Y.Kim
//    add koppenburg cut
//
// Revision 4.6  2010/05/03      T.Y.Kim
//    add koppenburg cut
//    add #2, #3, #4
//
// Revision 4.61  2010/06/07      T.Y.Kim
//    remove ptype_gg
//
// Revision 4.64  2010/06/29      T.Y.Kim
//    change ptype_gg -> ptype_pi0
//
// Revision 5.00  2010/09/07      T.Y.Kim
//    add ksfw
//
// Revision 5.01  2010/09/13      T.Y.Kim
//    add qr
//
// Revision 5.02  2010/09/14      T.Y.Kim
//    add dz
//
// Revision 4.71  2010/06/03      T.Y.Kim
//    (de,mbc) : (1.0, 5.19) -> (0.5, 5.20)  rec_b2gg
//
// Revision 4.72  2010/06/29      T.Y.Kim
//    change ptype_gg -> ptype_pi0
//
// Revision 4.73  2010/10/28      T.Y.Kim
//    add gmom
//
// Revision 4.74  2011/01/06      T.Y.Kim
//    change SCF, mc_b2gg=2  ->  mc_b2gg=2,3
//
// Revision 4.74  2011/01/12      T.Y.Kim
//    change SCF, mc_b2gg==2  ->  mc_b2gg=10,11
//    change true mc_b2gg==1  ->  mc_b2gg==1~9
//
// Revision 4.74  2011/02/24      T.Y.Kim
//    change SCF, mc_b2gg==2  ->  delete
//    change true mc_b2gg==1  ->  mc_b2gg==1~4
//
// Revision 4.75  2011/04/15      T.Y.Kim
//    tag for getting # of evt
//
// Revision 5.03  2011/04/15      T.Y.Kim
//    add ksfw, qr
//

  // e    11 
  // gam  22 
  // pi0 111 
  // eta 221 
  // B   511 


#include "b2gg.h"
#include "rec2mdst/ecltiming.h"
//
//
//
#if defined(BELLE_NAMESPACE)
namespace Belle {
#endif
//
//
//
void 
CLASS_NAME::event( BelleEvent* evptr, int* status )
{

  *status = 0;
  
  //bool n_evt_cut = false;  // tag for checking # of evt

  H[0]->accumulate( (float) 0.5 );
  // Belle event 
  if ( BsCouTab ( BELLE_EVENT ) ) {
    Belle_event_Manager& event_mgr = Belle_event_Manager::get_manager();
    Belle_event& event(event_mgr[0]);
    //
    EvtNo = event.EvtNo();
    Farm  = EvtNo>>28;
    Event = EvtNo-(Farm << 28);
  }
  // beam energy 
    ebeam = BeamEnergy::E_beam2();
    // ebeam = BeamEnergy::E_beam_corr();   20090724 correct by Y.Unno
  CMBoost = -BeamEnergy::CMBoost();

  // get ipProfile 
  m_ip = get_ip( 1, m_eip, m_eip_s );
  if( abs( m_ip.z() ) > 5.0 ) return;
  H[0]->accumulate( (float) 1.5 );
  //
  //
  //
  int mc_sig   = 0;
  if (ExpMC==2){
    for(int i=0; i<=1; i++){
      mc_genid[i] = -1;
    }
    mc_sig = mc_geninfo();
  }
 





  std::vector<Particle> list_gam;
  std::vector<Particle> list_gg;
  std::vector<Particle> list_b;

  list_gam.clear();
  list_gg.clear();
  list_b.clear();

  bool n_evt_cut = false;  // tag for checking # of evt
  //int n_evt_cut = 0;  // tag for checking # of evt

  //
  // prepare gamma list 
  //
  Mdst_gamma_Manager & gamma_mgr = Mdst_gamma_Manager::get_manager();
  for (std::vector<Mdst_gamma>::const_iterator iii = gamma_mgr.begin();
       iii != gamma_mgr.end(); iii++){
    //
    // extract timing information 
    //

    Mdst_ecl & m_ecl = iii->ecl();
    Mdst_ecl_Manager & m_ecl_m =  Mdst_ecl_Manager::get_manager();
    Mdst_ecl_aux_Manager & m_ecl_aux_m = Mdst_ecl_aux_Manager::get_manager();
    Mdst_ecl_aux_Manager::iterator m_ecl_aux_it = m_ecl_aux_m.begin();
    int tdccount = 0; 
    for (Mdst_ecl_Manager::iterator jjj = m_ecl_m.begin();
	 jjj != m_ecl_m.end(); jjj++){
      if (*jjj == m_ecl) {
	float prop2 = m_ecl_aux_it->property(2);
	int property2 = *reinterpret_cast<int*>(&prop2);
	float frc,frt;
	bool lcl,ltc;
	ecltiming::decode(property2,tdccount,frc,frt,lcl,ltc);
      }
      m_ecl_aux_it++;
    }
    //
    // e9/e25 
    //
    Mdst_ecl_aux &ecl_aux = m_ecl_aux_m(Panther_ID(m_ecl.get_ID()));
    //
    // put timing information to particle 
    //
    Particle ppp(*iii);
    //    ppp.userInfo(*(new UserInfo()));
    //    UserInfo & info = dynamic_cast<UserInfo&>(ppp.userInfo());
    //    info.timing(tdccount);
    setUserInfo(ppp);
    dynamic_cast<UserInfo&>(ppp.userInfo()).timing(tdccount);
    dynamic_cast<UserInfo&>(ppp.userInfo()).e9oe25(ecl_aux.e9oe25());
    list_gam.push_back( ppp );
    //
    //
    //
  }
  if (list_gam.size()<2){ return; }
  H[0]->accumulate( (float) 2.5 );






  //
  //  pi0,eta by Mdst_gamma 
  //
  combination(list_gg, ptype_pi0, list_gam, list_gam);   //ptype_gg -> ptype_pi0

  for( std::vector<Particle>::iterator iii = list_gg.begin();
       iii != list_gg.end(); ++iii ){
 
    Particle ptcl_gg(*iii);
    Particle ptcl_gam01;
    Particle ptcl_gam02;
    if (iii->child(0).p().vect().mag()>iii->child(1).p().vect().mag()){
      ptcl_gam01 = iii->child(0);
      ptcl_gam02 = iii->child(1);
    }else{
      ptcl_gam01 = iii->child(1);
      ptcl_gam02 = iii->child(0);
    }

    const Mdst_gamma gam01 = ptcl_gam01.mdstGamma();
    const Mdst_gamma gam02 = ptcl_gam02.mdstGamma();

    HepLorentzVector p4lab_gam01 = ptcl_gam01.p();
    HepLorentzVector p4lab_gam02 = ptcl_gam02.p();
    HepLorentzVector p4lab_gam12 = p4lab_gam01 + p4lab_gam02;
    double mass_g12 = p4lab_gam12.mag();

    double gam01_e = p4lab_gam01.e();
    double gam02_e = p4lab_gam02.e();
    double gg_lang = cos(p4lab_gam02.angle(p4lab_gam01.vect()));
    float gam01_theta = p4lab_gam01.theta()*180.0/CC_PI;
    float gam02_theta = p4lab_gam02.theta()*180.0/CC_PI;
 

    HepLorentzVector p4lab_gg = ptcl_gg.p();
    HepLorentzVector p4cms_gg = p4lab_gg;
    p4cms_gg.boost(CMBoost);

    double gg_mass = ptcl_gg.mass();
    float gg_le   = p4lab_gg.e();
    float gg_ce   = p4cms_gg.e();
    float gg_lmom = p4lab_gg.vect().mag();
    float gg_cmom = p4cms_gg.vect().mag();
    float gg_lcos = p4lab_gg.vect().cosTheta();
    float gg_ccos = p4cms_gg.vect().cosTheta();
    //float gg_chi2 = ptcl_gg.chisq();
 

  
    //Barrel, endcap e, theta cut
    int flag_ecut = 0; 

    if ( gam01_theta< 32 || gam01_theta>130 ) {
      if ( gam01_e<0.100 ) { flag_ecut=1;}
    }else{
      if ( gam01_e<0.050 ) { flag_ecut=1;}
    }
    if ( gam02_theta<32 || gam02_theta>130) {
      if ( gam02_e<0.100 ) { flag_ecut=1;}
    }else{
      if ( gam02_e<0.050 ) { flag_ecut=1;}
    }

    if (flag_ecut==0){

    // gg_flag == 2   pi0
    // gg_flag == 4   eta
    // gg_flag == 0   others

    int gg_flag01 = pi0etaveto ( gam01 );
    int gg_flag02 = pi0etaveto ( gam02 );
    int gg_flag = gg_flag01 + gg_flag02;

    Mdst_gamma g1_pi0;
    Mdst_gamma g2_pi0;
    Mdst_gamma g4_pi0;
    Mdst_gamma g5_pi0;
    Mdst_gamma g1_eta;
    Mdst_gamma g2_eta;
    Mdst_gamma g4_eta;
    Mdst_gamma g5_eta;
    
    double m_pi1, m_pi2, m_pi4, m_pi5, m_et1, m_et2, m_et4, m_et5;

    double prob_pi1 = Pi0_Prob(gg_mass, gam01_e, gam01_theta);    
    double prob_pi2 = Pi0_Prob(gg_mass, gam02_e, gam02_theta);
    double prob_et1 = Eta_Prob(gg_mass, gam01_e, gam01_theta);
    double prob_et2 = Eta_Prob(gg_mass, gam02_e, gam02_theta);    

    double prob_pi0 = prob_pi1 + prob_pi2;
    double prob_eta = prob_et1 + prob_et2;

    double clos_pi1 = Closest_Pi0_Probability(gam01, g1_pi0, m_pi1, gam01);
    double clos_pi2 = Closest_Pi0_Probability(gam02, g2_pi0, m_pi2, gam02);
    double clos_et1 = Closest_Eta_Probability(gam01, g1_eta, m_et1, gam01);
    double clos_et2 = Closest_Eta_Probability(gam02, g2_eta, m_et2, gam02);

    double clos_pi4 = Closest_Pi0_Probability(gam01, g4_pi0, m_pi4, gam02);
    double clos_pi5 = Closest_Pi0_Probability(gam02, g5_pi0, m_pi5, gam01);
    double clos_et4 = Closest_Eta_Probability(gam01, g4_eta, m_et4, gam02);
    double clos_et5 = Closest_Eta_Probability(gam02, g5_eta, m_et5, gam01);

    double clos_pi0 = clos_pi1 + clos_pi2;
    double clos_eta = clos_et1 + clos_et2;
    double clos_pi3 = clos_pi4 + clos_pi5;
    double clos_et3 = clos_et4 + clos_et5;

    double m_pi0 = m_pi1 + m_pi2;
    double m_eta = m_et1 + m_et2;
    double m_pi3 = m_pi4 + m_pi5;
    double m_et3 = m_et4 + m_et5;

    //
    //  MC
    //
    int mc_idhep_gg = 0;
    if(ExpMC==2){
      mc_idhep_gg = mc_gg(ptcl_gg);
    }

    T[2]->column( "gg_mass", (float) gg_mass );
    T[2]->column( "gg_le  ", (float) gg_le );
    T[2]->column( "gg_ce  ", (float) gg_ce );
    T[2]->column( "gg_lmom", (float) gg_lmom );//mass constraint
    T[2]->column( "gg_cmom", (float) gg_cmom );//mass constraint
    T[2]->column( "gg_lcos", (float) gg_lcos );//mass constraint
    T[2]->column( "gg_ccos", (float) gg_ccos );//mass constraint
    //T[2]->column( "pi0_chi2", (float) pi0_chi2 );   

    T[2]->column( "veto1",   (float) gg_flag01);
    T[2]->column( "veto2",   (float) gg_flag02); 
    T[2]->column( "veto",    (float) gg_flag);

    T[2]->column( "g1_pi0",   (float) g1_pi0 );
    T[2]->column( "g1_eta",   (float) g1_eta );  
    T[2]->column( "g2_pi0",   (float) g2_pi0 );
    T[2]->column( "g2_eta",   (float) g2_eta ); 
    T[2]->column( "g4_pi0",   (float) g4_pi0 );
    T[2]->column( "g4_eta",   (float) g4_eta ); 
    T[2]->column( "g5_pi0",   (float) g5_pi0 );
    T[2]->column( "g5_eta",   (float) g5_eta ); 
      
    T[2]->column( "m_pi1",    (float) m_pi1 );
    T[2]->column( "m_pi2",    (float) m_pi2 );
    T[2]->column( "m_pi4",    (float) m_pi4 );
    T[2]->column( "m_pi5",    (float) m_pi5 );
    T[2]->column( "m_et1",    (float) m_et1 ); 
    T[2]->column( "m_et2",    (float) m_et2 ); 
    T[2]->column( "m_et4",    (float) m_et4 ); 
    T[2]->column( "m_et5",    (float) m_et5 ); 
    T[2]->column( "m_pi0",    (float) m_pi0 );
    T[2]->column( "m_eta",    (float) m_eta ); 
    T[2]->column( "m_pi3",    (float) m_pi3 );
    T[2]->column( "m_et3",    (float) m_et3 ); 

    T[2]->column( "prob_pi0", (float) prob_pi0 );
    T[2]->column( "prob_eta", (float) prob_eta );
    T[2]->column( "prob_pi1", (float) prob_pi1 );
    T[2]->column( "prob_pi2", (float) prob_pi2 );
    T[2]->column( "prob_et1", (float) prob_et1 );
    T[2]->column( "prob_et2", (float) prob_et2 );
      
    T[2]->column( "clos_pi0", (float) clos_pi0 );
    T[2]->column( "clos_eta", (float) clos_eta );
    T[2]->column( "clos_pi1", (float) clos_pi1 );
    T[2]->column( "clos_pi2", (float) clos_pi2 );
    T[2]->column( "clos_pi3", (float) clos_pi3 );
    T[2]->column( "clos_pi4", (float) clos_pi4 );
    T[2]->column( "clos_pi5", (float) clos_pi5 );
    T[2]->column( "clos_et1", (float) clos_et1 );
    T[2]->column( "clos_et2", (float) clos_et2 );
    T[2]->column( "clos_et3", (float) clos_et3 );
    T[2]->column( "clos_et4", (float) clos_et4 );
    T[2]->column( "clos_et5", (float) clos_et5 );
      
    T[2]->column( "mc_gg",   (float) mc_idhep_gg );
    T[2]->dumpData();

    }

  }






  combination(list_b, ptype_b, list_gam, list_gam);

  int n_cand = rec_b2gg(list_b);
  
  if (n_cand==0){ return; }
  T[0]->column(    "ncand",      (float) n_cand );
  H[0]->accumulate( (float) 3.5 );







  //
  // Bhabha veto trigger information 
  //
  float t01 = 0;
  float t02 = 0;
  float t0n = 0;
  std::vector<float> t0 = TriggerTiming(TrgBit::csi_bb);
  t0n = 0;
  for( std::vector<float>::iterator iii = t0.begin();
       iii != t0.end(); ++iii ){
    t0n++;
    if      (t0n==1){ t01 = *iii; }
    else if (t0n==2){ t02 = *iii; }
  }
    T[1]->column( "expno",   (float) ExpNo );
    T[1]->column( "runno",   (float) RunNo );
    T[1]->column( "evtno",   (float) EvtNo );
    T[1]->column( "t01",     (float) t01 );
    T[1]->column( "t02",     (float) t02 );
    T[1]->column( "t0n",     (float) t0n );
    T[1]->dumpData();


 


    //  combination(list_gg, ptype_gg, list_gam, list_gam);


  for( std::vector<Particle>::iterator iii = list_b.begin();
       iii != list_b.end(); ++iii ){

    Particle ptcl_b(*iii);
    Particle ptcl_gam01;
    Particle ptcl_gam02;
    if (iii->child(0).p().vect().mag()>iii->child(1).p().vect().mag()){
      ptcl_gam01 = iii->child(0);
      ptcl_gam02 = iii->child(1);
    }else{
      ptcl_gam01 = iii->child(1);
      ptcl_gam02 = iii->child(0);
    }

    HepLorentzVector p4lab_b = ptcl_b.p();
    HepLorentzVector p4cms_b = p4lab_b;
    p4cms_b.boost(CMBoost);

    //
    // vertex 
    //
    double vz_sig = vertex_sig(*iii);
    double vz_tag = vertex_tag(*iii);
    double dz = vz_sig - vz_tag;
    //
    // fbtag 
    //
    double qr = fbtg(*iii);
    //
    //
    //

    double mbc2 = pow( ebeam, 2. ) - p4cms_b.vect().mag2();
    if (mbc2<0){continue;}
    double mbc  = sqrt( mbc2 );
    double de   = p4cms_b.t() - ebeam;
    if ( abs(de)>2.0 ||mbc<5.20 ){ continue; }
    //
    //
    //
    const Mdst_gamma gam01 = ptcl_gam01.mdstGamma();
    const Mdst_gamma gam02 = ptcl_gam02.mdstGamma();

    Hep3Vector p3lab_gam01(gam01.px(), gam01.py(), gam01.pz());
    Hep3Vector p3lab_gam02(gam02.px(), gam02.py(), gam02.pz());
    HepLorentzVector p4lab_gam01(p3lab_gam01,p3lab_gam01.mag());
    HepLorentzVector p4lab_gam02(p3lab_gam02,p3lab_gam02.mag());

    HepLorentzVector p4cms_gam01 = p4lab_gam01;
    HepLorentzVector p4cms_gam02 = p4lab_gam02;
    p4cms_gam01.boost(CMBoost);
    p4cms_gam02.boost(CMBoost);    

    HepLorentzVector p4lab_gam12 = p4lab_gam01 + p4lab_gam02;
    double gg_mass = p4lab_gam12.mag();
  
    float gam01_e = p4lab_gam01.e();
    float gam02_e = p4lab_gam02.e();

    float gam01_theta = p4lab_gam01.theta()*180.0/CC_PI;
    float gam02_theta = p4lab_gam02.theta()*180.0/CC_PI;
    double gg_lphi = CC_PI - fabs(p4lab_gam01.phi()-p4lab_gam02.phi());
    double gg_cphi = CC_PI - fabs(p4cms_gam01.phi()-p4cms_gam02.phi());
    double gg_lang = cos(p4lab_gam02.angle(p4lab_gam01.vect()));
    double gg_cang = cos(p4cms_gam02.angle(p4cms_gam01.vect()));
 
    //Barrel, endcap e, theta cut
    int flag_ecut = 0; 

    if ( gam01_theta< 32 || gam01_theta>130 ) {
      if ( gam01_e<0.100 ) { flag_ecut=1;}
    }else{
      if ( gam01_e<0.050 ) { flag_ecut=1;}
    }
    if ( gam02_theta<32 || gam02_theta>130) {
      if ( gam02_e<0.100 ) { flag_ecut=1;}
    }else{
      if ( gam02_e<0.050 ) { flag_ecut=1;}
    }

    if (flag_ecut==0){

    H[10]->accumulate( (float) gg_mass );


    // gg_flag == 2   pi0
    // gg_flag == 4   eta
    // gg_flag == 0   others

    int gg_flag01 = pi0etaveto ( gam01 );
    int gg_flag02 = pi0etaveto ( gam02 );
    int gg_flag = gg_flag01 + gg_flag02;
    

    Mdst_gamma g1_pi0;
    Mdst_gamma g2_pi0;
    Mdst_gamma g4_pi0;
    Mdst_gamma g5_pi0;
    Mdst_gamma g1_eta;
    Mdst_gamma g2_eta;
    Mdst_gamma g4_eta;
    Mdst_gamma g5_eta;
    
    double m_pi1, m_pi2, m_pi4, m_pi5, m_et1, m_et2, m_et4, m_et5;

    double prob_pi1 = Pi0_Prob(gg_mass, gam01_e, gam01_theta);    
    double prob_pi2 = Pi0_Prob(gg_mass, gam02_e, gam02_theta);
    double prob_et1 = Eta_Prob(gg_mass, gam01_e, gam01_theta);
    double prob_et2 = Eta_Prob(gg_mass, gam02_e, gam02_theta);    

    double prob_pi0 = prob_pi1 + prob_pi2;
    double prob_eta = prob_et1 + prob_et2;

    double clos_pi1 = Closest_Pi0_Probability(gam01, g1_pi0, m_pi1, gam01);
    double clos_pi2 = Closest_Pi0_Probability(gam02, g2_pi0, m_pi2, gam02);
    double clos_et1 = Closest_Eta_Probability(gam01, g1_eta, m_et1, gam01);
    double clos_et2 = Closest_Eta_Probability(gam02, g2_eta, m_et2, gam02);

    double clos_pi4 = Closest_Pi0_Probability(gam01, g4_pi0, m_pi4, gam02);
    double clos_pi5 = Closest_Pi0_Probability(gam02, g5_pi0, m_pi5, gam01);
    double clos_et4 = Closest_Eta_Probability(gam01, g4_eta, m_et4, gam02);
    double clos_et5 = Closest_Eta_Probability(gam02, g5_eta, m_et5, gam01);

    double clos_pi0 = clos_pi1 + clos_pi2;
    double clos_eta = clos_et1 + clos_et2;
    double clos_pi3 = clos_pi4 + clos_pi5;
    double clos_et3 = clos_et4 + clos_et5;

    double m_pi0 = m_pi1 + m_pi2;
    double m_eta = m_et1 + m_et2;
    double m_pi3 = m_pi4 + m_pi5;
    double m_et3 = m_et4 + m_et5;

    //
    // ksfw moments
    //
    ksfwmoments km( *iii, ebeam, CMBoost);

    //
    // MC 
    //
    int mc_b2gg = 0;
    int mc_g1 = 0;
    int mc_g2 = 0;
    int mc_g1_mom = 0;
    int mc_g2_mom = 0;
    int mc_g1_gmom = 0;
    int mc_g2_gmom = 0;



    for( std::vector<Particle>::iterator jjj = list_gg.begin();
	 jjj != list_gg.end(); ++jjj ){
      
      Particle ptcl_gg(*jjj);
      
      int mc_idhep_gg = 0;
      if(ExpMC==2){
	mc_idhep_gg = mc_gg(ptcl_gg);
      }
      
      T[0]->column( "mc_gg",    (float) mc_idhep_gg );
    }
    

    
    if(ExpMC==2){

      mc_g1 = mc_idhep(ptcl_gam01);
      mc_g2 = mc_idhep(ptcl_gam02);

      mc_b2gg   = u_mc_type(ptcl_gam01,ptcl_gam02);
      // mother 
      mc_g1_mom = mc_idhep_mother(&ptcl_gam01.mdstGamma());
      mc_g2_mom = mc_idhep_mother(&ptcl_gam02.mdstGamma());
      mc_g1_gmom = mc_idhep_gmother(&ptcl_gam01.mdstGamma());
      mc_g2_gmom = mc_idhep_gmother(&ptcl_gam02.mdstGamma());

    }
    //
    // timing 
    //
    int timing_g01 = dynamic_cast<UserInfo&>(ptcl_gam01.userInfo()).timing();
    int timing_g02 = dynamic_cast<UserInfo&>(ptcl_gam02.userInfo()).timing();
    //
    // ntuple 
    //
    T[0]->column( "expno",   (float) ExpNo );
    T[0]->column( "runno",   (float) RunNo );
    T[0]->column( "evtno",   (float) EvtNo );
    T[0]->column( "t01",     (float)   t01 );
    T[0]->column( "t02",     (float)   t02 );
    T[0]->column( "t0n",     (float)   t0n );
    T[0]->column( "ebeam",   (float) ebeam );
    //
    T[0]->column( "mbc",     (float) mbc );
    T[0]->column( "de",      (float) de );

    T[0]->column( "g1_lmom", (float) p4lab_gam01.vect().mag() );
    T[0]->column( "g2_lmom", (float) p4lab_gam02.vect().mag() );
    T[0]->column( "g1_cmom", (float) p4cms_gam01.vect().mag() );
    T[0]->column( "g2_cmom", (float) p4cms_gam02.vect().mag() );
    T[0]->column( "g1_lcos", (float) p4lab_gam01.vect().cosTheta() );
    T[0]->column( "g2_lcos", (float) p4lab_gam02.vect().cosTheta() );
    T[0]->column( "g1_ccos", (float) p4cms_gam01.vect().cosTheta() );
    T[0]->column( "g2_ccos", (float) p4cms_gam02.vect().cosTheta() );
    T[0]->column( "g1_e925", (float) dynamic_cast<UserInfo&>(ptcl_gam01.userInfo()).e9oe25() );
    T[0]->column( "g2_e925", (float) dynamic_cast<UserInfo&>(ptcl_gam02.userInfo()).e9oe25() );
    T[0]->column( "gg_lang", (float) gg_lang );
    T[0]->column( "gg_cang", (float) gg_cang );
    T[0]->column( "gg_lphi", (float) gg_lphi );
    T[0]->column( "gg_cphi", (float) gg_cphi );
    T[0]->column( "gg_mass",  (float) gg_mass );

    T[0]->column( "gam01_e", (float) gam01_e );
    T[0]->column( "gam02_e", (float) gam02_e );

    T[0]->column( "g1_t",     (int)   timing_g01 );
    T[0]->column( "g2_t",     (int)   timing_g02 );

    T[0]->column( "mc_b2gg",  (float) mc_b2gg );
    T[0]->column( "mc_g1",    (float) mc_g1 );
    T[0]->column( "mc_g2",    (float) mc_g2 );
    T[0]->column( "mc_g1m",   (float) mc_g1_mom );
    T[0]->column( "mc_g2m",   (float) mc_g2_mom );
    T[0]->column( "mc_g1gm",   (float) mc_g1_gmom );
    T[0]->column( "mc_g2gm",   (float) mc_g2_gmom );

    T[0]->column( "g1_pi0",   (float) g1_pi0 );
    T[0]->column( "g1_eta",   (float) g1_eta );  
    T[0]->column( "g2_pi0",   (float) g2_pi0 );
    T[0]->column( "g2_eta",   (float) g2_eta ); 
    T[0]->column( "g4_pi0",   (float) g4_pi0 );
    T[0]->column( "g4_eta",   (float) g4_eta ); 
    T[0]->column( "g5_pi0",   (float) g5_pi0 );
    T[0]->column( "g5_eta",   (float) g5_eta ); 

    T[0]->column( "veto1",   (float) gg_flag01);
    T[0]->column( "veto2",   (float) gg_flag02); 
    T[0]->column( "veto",    (float) gg_flag);

    T[0]->column( "m_pi1",    (float) m_pi1 );
    T[0]->column( "m_pi2",    (float) m_pi2 );
    T[0]->column( "m_pi4",    (float) m_pi4 );
    T[0]->column( "m_pi5",    (float) m_pi5 );
    T[0]->column( "m_et1",    (float) m_et1 ); 
    T[0]->column( "m_et2",    (float) m_et2 ); 
    T[0]->column( "m_et4",    (float) m_et4 ); 
    T[0]->column( "m_et5",    (float) m_et5 ); 
    T[0]->column( "m_pi0",    (float) m_pi0 );
    T[0]->column( "m_eta",    (float) m_eta ); 
    T[0]->column( "m_pi3",    (float) m_pi3 );
    T[0]->column( "m_et3",    (float) m_et3 ); 

    T[0]->column( "prob_pi0", (float) prob_pi0 );
    T[0]->column( "prob_eta", (float) prob_eta );
    T[0]->column( "prob_pi1", (float) prob_pi1 );
    T[0]->column( "prob_pi2", (float) prob_pi2 );
    T[0]->column( "prob_et1", (float) prob_et1 );
    T[0]->column( "prob_et2", (float) prob_et2 );

    T[0]->column( "clos_pi0", (float) clos_pi0 );
    T[0]->column( "clos_eta", (float) clos_eta );
    T[0]->column( "clos_pi1", (float) clos_pi1 );
    T[0]->column( "clos_pi2", (float) clos_pi2 );
    T[0]->column( "clos_pi3", (float) clos_pi3 );
    T[0]->column( "clos_pi4", (float) clos_pi4 );
    T[0]->column( "clos_pi5", (float) clos_pi5 );
    T[0]->column( "clos_et1", (float) clos_et1 );
    T[0]->column( "clos_et2", (float) clos_et2 );
    T[0]->column( "clos_et3", (float) clos_et3 );
    T[0]->column( "clos_et4", (float) clos_et4 );
    T[0]->column( "clos_et5", (float) clos_et5 );

    //T[0]->column( "lundgmo1", (float) gam01.genHepevt().mother() ? gam01.genHepevt().mother().idhep() : 0.);
    //T[0]->column( "lundgmo2", (float) gam02.genHepevt().mother() ? gam02.genHepevt().mother().idhep() : 0.);

    T[0]->column( "cosb",     (float) p4cms_b.vect().cosTheta() );
    T[0]->column( "dz",       (float) dz );
    T[0]->column( "qr",       (float) qr );

    km.usefinal(0);
    T[0]->column("k0mm2",   km.mm2());
    T[0]->column("k0et",    km.et());
    T[0]->column("k0hso00", km.Hso(0, 0));
    T[0]->column("k0hso01", km.Hso(0, 1));
    T[0]->column("k0hso02", km.Hso(0, 2));
    T[0]->column("k0hso03", km.Hso(0, 3));
    T[0]->column("k0hso04", km.Hso(0, 4));
    T[0]->column("k0hso10", km.Hso(1, 0));
    T[0]->column("k0hso12", km.Hso(1, 2));
    T[0]->column("k0hso14", km.Hso(1, 4));
    T[0]->column("k0hso20", km.Hso(2, 0));
    T[0]->column("k0hso22", km.Hso(2, 2));
    T[0]->column("k0hso24", km.Hso(2, 4));
    T[0]->column("k0hoo0",  km.Hoo(0));
    T[0]->column("k0hoo1",  km.Hoo(1));
    T[0]->column("k0hoo2",  km.Hoo(2));
    T[0]->column("k0hoo3",  km.Hoo(3));
    T[0]->column("k0hoo4",  km.Hoo(4));

    km.usefinal(1);
    T[0]->column("k1mm2",   km.mm2());
    T[0]->column("k1et",    km.et()); 
    T[0]->column("k1hso00", km.Hso(0, 0));
    T[0]->column("k1hso01", km.Hso(0, 1));
    T[0]->column("k1hso02", km.Hso(0, 2));
    T[0]->column("k1hso03", km.Hso(0, 3));
    T[0]->column("k1hso04", km.Hso(0, 4));
    T[0]->column("k1hso10", km.Hso(1, 0));
    T[0]->column("k1hso12", km.Hso(1, 2));
    T[0]->column("k1hso14", km.Hso(1, 4));
    T[0]->column("k1hso20", km.Hso(2, 0));
    T[0]->column("k1hso22", km.Hso(2, 2));
    T[0]->column("k1hso24", km.Hso(2, 4));
    T[0]->column("k1hoo0",  km.Hoo(0));
    T[0]->column("k1hoo1",  km.Hoo(1));
    T[0]->column("k1hoo2",  km.Hoo(2));
    T[0]->column("k1hoo3",  km.Hoo(3));
    T[0]->column("k1hoo4",  km.Hoo(4));

    T[0]->dumpData();

    //if ( abs(de)<0.5 && mbc>5.20 ){ n_evt_cut = true ;}  // tag for checking # of evt


    } // flag_ecut
    
    if ( abs(de)<0.5 && mbc>5.20 ){ n_evt_cut = true ;}  // tag for checking # of evt
    //if ( abs(de)<0.5 && mbc>5.20 ){ n_evt_cut = 1 ;}  // tag for checking # of evt

  } //loop for list_b

    if( n_evt_cut == true ) {  H[11]->accumulate( (float) 4.5 ); } // check # of evt




  //
  //  pi0,eta veto by Mdst_gamma
  //
  for( std::vector<Particle>::iterator iii = list_b.begin();
       iii != list_b.end(); ++iii ){
 
    Particle ptcl_b(*iii);
 
    HepLorentzVector p4lab_b = ptcl_b.p();
    HepLorentzVector p4cms_b = p4lab_b;
    p4cms_b.boost(CMBoost);

    double mbc2 = pow( ebeam, 2. ) - p4cms_b.vect().mag2();
    if (mbc2<0){continue;}
    double mbc  = sqrt( mbc2 );
    double de   = p4cms_b.t() - ebeam;
    if ( abs(de)>2.0 ||mbc<5.20 ){ continue; }
    //
    //
    //
    Particle ptcl_gam01;
    Particle ptcl_gam02;
    if (iii->child(0).p().vect().mag()>iii->child(1).p().vect().mag()){
      ptcl_gam01 = iii->child(0);
      ptcl_gam02 = iii->child(1);
    }else{
      ptcl_gam01 = iii->child(1);
      ptcl_gam02 = iii->child(0);
    }

    const Mdst_gamma gam01 = ptcl_gam01.mdstGamma();
    const Mdst_gamma gam02 = ptcl_gam02.mdstGamma();

    Hep3Vector p3lab_gam01(gam01.px(), gam01.py(), gam01.pz());
    Hep3Vector p3lab_gam02(gam02.px(), gam02.py(), gam02.pz());
    HepLorentzVector p4lab_gam01(p3lab_gam01,p3lab_gam01.mag());
    HepLorentzVector p4lab_gam02(p3lab_gam02,p3lab_gam02.mag());
  
    float gam01_e = p4lab_gam01.e();
    float gam02_e = p4lab_gam02.e();
    double gg_lang = cos(p4lab_gam02.angle(p4lab_gam01.vect()));
    float gam01_theta = p4lab_gam01.theta()*180.0/CC_PI;
    float gam02_theta = p4lab_gam02.theta()*180.0/CC_PI;
 

    //Barrel, endcap e, theta cut
    int flag_ecut = 0; 

    if ( gam01_theta< 32 || gam01_theta>130 ) {
      if ( gam01_e<0.100 ) { flag_ecut=1;}
    }else{
      if ( gam01_e<0.050 ) { flag_ecut=1;}
    }
    if ( gam02_theta<32 || gam02_theta>130) {
      if ( gam02_e<0.100 ) { flag_ecut=1;}
    }else{
      if ( gam02_e<0.050 ) { flag_ecut=1;}
    }

    if (flag_ecut==0){

    // gg_flag == 2   pi0
    // gg_flag == 4   eta
    // gg_flag == 0   others

    int gg_flag01 = pi0etaveto ( gam01 );
    int gg_flag02 = pi0etaveto ( gam02 );
    int gg_flag = gg_flag01 + gg_flag02;
    
    //
    //
    //
    T[3]->column( "expno",   (float) ExpNo );
    T[3]->column( "runno",   (float) RunNo );
    T[3]->column( "evtno",   (float) EvtNo );
    T[3]->column( "mbc",     (float) mbc );
    T[3]->column( "de",      (float) de );

    T[3]->column( "veto1",   (float) gg_flag01);
    T[3]->column( "veto2",   (float) gg_flag02); 
    T[3]->column( "veto",    (float) gg_flag);

    T[3]->dumpData();
    }
  }  //pi0etaveto loop
    
} //event






 
double CLASS_NAME::fbtg(Particle &p){

 
  /* example for baseline method
  Particle Bcp; // reconstruct Bcp using Particle class
  Hamlet hamlet; // define Hamlet object
  hamlet.setBcp( Bcp ); // set Bcp particle (and fill Brecon table)
  // get result
  int flavor = hamlet.flavor(); // +1 : for tag-side being B0 or B+
                                // -1 : for tag-side being B0-bar or B-
                                //  0 : does not specify a b-flavor
                                // -2 : does not specify a b-flavor
				//      (no leptons and kaons are found) 
  int imode  = hamlet.imode(); // 1 : tag a b-flavor with high-electron method
                               // 2 : tag a b-flavor with high-muon method
                               // 4 : tag a b-flavor with kaon method
                               // others : does not specify a b-flavor
  */

  /* example for MDLH
  Hamlet hamlet; // Define Hamlet
  Particle Bcp=...; // Reconstruct CP side B meson
  hamlet.setBcp(Bcp); // Set Bcp
  hamlet.setTagMethod(Hamlet::MULT_DIM_LH0); // Set hamlet to MDLH method
  double qr = hamlet.q(); // get q*r
  // get the reference of MDLH class (see header file for its member function)
  const Fbtag_MultDimLikelihood0 & mdlh0(hamlet.fbtg_mult_dim_likelihood0());
  */
  
  Hamlet hamlet;
  hamlet.setBcp(p);
  //  hamlet.tag( Hamlet::BASELINE );
  // hamlet.tag( (Hamlet::tag_method)19 );
  hamlet.tag((Hamlet::tag_method)18);
  /*
  @param MULT_DIM_LH   (18) tag by mult-dimensional likelihood method
  This was (19) up to b20020424_1007.
  */
  double qr= hamlet.q();
  /*
  T->column( "qr", (float)hamlet.q() );
  T->column( "flavor", (float)hamlet.flavor() );
  T->column( "imode", (float)hamlet.imode() );
  T->column( "wtag", (float)hamlet.wtag() );
  */
  return qr;
    
}









double 
CLASS_NAME::vertex_sig( Particle &ppp)
{
  //
  kvertexfitter kv;
  //
  int n_used_tracks = 0;
  for( int i=0; i<ppp.relation().nFinalStateParticles(); i++ ){
    const Particle &p(ppp.relation().finalStateParticle(i));
    // charged track? 
    if( !p.mdstCharged() ) continue;
    if( !enough_svd_hit(p, 1, 2) ) continue;
    // ok.  add to kfitter's particle bank
    addTrack2fit( kv, p );
    n_used_tracks++;
  }
  //
  if( IpProfile::b_life_smeared_usable() ){
    addBeam2fit( kv, m_ip, m_eip_s );
  }
  //
  //
  //
  const int error = kv.fit();
  double chisq  = 1e9;
  double ndf    = 999.0;
  double cl     = -1.0;
  double chisq2 = 1e9;
  double ndf2   = 999.0;
  double cl2    = -1.0;
  if( !error ){
    ppp.momentum().decayVertex( kv.vertex(), kv.errVertex() );
    chisq  = kv.chisq();
    ndf    = kv.dgf();
    cl     = kv.cl();
    chisq2 = kv.chisq_woip();
    ndf2   = kv.dgf_woip();
    cl2    = kv.cl_woip();
  } else {
    ppp.momentum().decayVertex( HepPoint3D( 999., 999., 999. ) );
  }
  T[0]->column( "cvx",    (float)ppp.momentum().decayVertex().x() );
  T[0]->column( "cvy",    (float)ppp.momentum().decayVertex().y() );
  T[0]->column( "cvz",    (float)ppp.momentum().decayVertex().z() );
  T[0]->column( "cevz",   (float)sqrt( ppp.momentum().dDecayVertex().fast(3,3) ) );
  T[0]->column( "cchisq", (float)chisq );
  T[0]->column( "cndf",   (float)ndf );
  T[0]->column( "ccl",    (float)cl );
  T[0]->column( "cchisq2",(float)chisq2 );
  T[0]->column( "cndf2",  (float)ndf2 );
  T[0]->column( "ccl2",   (float)cl2 );
  T[0]->column( "cntrk",  (float)n_used_tracks);

  return ppp.momentum().decayVertex().z();
}








double 
CLASS_NAME::vertex_tag( const Particle &ppp )
{
  Particle p(ppp);
  //
  std::vector<Particle> tag_list;
  mktaglist( p, tag_list);
  //
  TagVK v_tag;

  v_tag.setdefault( p, ppp.momentum().decayVertex());
  //  v_tag.setdefault( p, p.momentum().decayVertex());
  //
  int remained( 0 );
  for( std::vector<Particle>::iterator it = tag_list.begin();
      it != tag_list.end(); ++it ){
    v_tag.push_back( &(*it) );
    remained++;
  }
  //
  HepPoint3D vtx(-999.0, -999.0, -999.0);
  HepSymMatrix evtx( 3, 0 );
  double chisq  = -1.0;
  double ndf    = 999.0;
  double cl     = -1.0;
  double chisq2 = -1.0;
  double ndf2   = 999.0;
  double cl2    = -1;

  if( !v_tag.fit() ){
    vtx    = v_tag.vtx();
    evtx   = v_tag.errVtx();
    chisq  = v_tag.chisq();
    chisq2 = v_tag.chisq_woip();
    ndf    = v_tag.ndf();
    ndf2   = v_tag.ndf_woip();
    cl     = v_tag.cl();
    cl2    = v_tag.cl_woip();
  }
  T[0]->column( "tvx",     (float)vtx.x() );
  T[0]->column( "tvy",     (float)vtx.y() );
  T[0]->column( "tvz",     (float)vtx.z() );
  T[0]->column( "tevz",    (float)sqrt( evtx.fast(3,3) ) );
  T[0]->column( "tchisq",  (float)chisq );
  T[0]->column( "tndf",    (float)ndf );
  T[0]->column( "tcl",     (float)cl );
  T[0]->column( "tchisq2", (float)chisq2 );
  T[0]->column( "tndf2",   (float)ndf2 );
  T[0]->column( "tcl2",    (float)cl2 );

  return vtx.z();
}










void 
CLASS_NAME::mktaglist( const Particle &p, 
		       std::vector<Particle> &taglist)
{
  std::vector<Particle> k_p, k_m, pi_p, pi_m;
  makeKPi( k_p, k_m, pi_p, pi_m, 0 );

  taglist = pi_p;
  taglist.insert( taglist.end(), pi_m.begin(), pi_m.end() );

  for( int i=0; i < p.relation().nFinalStateParticles(); ++i ){
    const Particle &child(p.relation().finalStateParticle(i));
    if( !child.mdstCharged() ) continue;
    removeParticle( taglist, p.relation().finalStateParticle(i));
  }
}









int CLASS_NAME::pi0etaveto(const Mdst_gamma &gam01)
{
  Mdst_gamma_Manager & gamma_mgr = Mdst_gamma_Manager::get_manager();

  int gg_flag = 0;
 
  for (std::vector<Mdst_gamma>::iterator kkk = gamma_mgr.begin();
       kkk != gamma_mgr.end(); kkk++){
    
    const Mdst_gamma &gam02 = *kkk;
    

    if(gam01.get_ID() == gam02.get_ID()) continue;    
    
    //    Particle ptcl_gam01(gam01, Ptype("GAMM"));
    //    Particle ptcl_gam02(gam02, Ptype("GAMM"));
    Particle ptcl_gam01(gam01);
    Particle ptcl_gam02(gam02);

    HepLorentzVector p4lab_gam01 = ptcl_gam01.p();
    HepLorentzVector p4lab_gam02 = ptcl_gam02.p();
    HepLorentzVector p4lab_gg    = p4lab_gam01 + p4lab_gam02;
    float gg_mass = p4lab_gg.mag();

    float gam01_e = p4lab_gam01.e();
    float gam02_e = p4lab_gam02.e();

    //H[10]->accumulate( gg_mass );
    //H[11]->accumulate( gg_mass );    
    
    //pi0 veto
    if ( (fabs(gg_mass-0.135) < 0.018) && (gam01_e>0.03 && gam02_e > 0.03) ) { gg_flag = 1;}
    //eta veto  
    else if ( (fabs(gg_mass-0.548) < 0.032 ) && (gam01_e>0.20 && gam02_e > 0.20)) { gg_flag = 2;}

    else { gg_flag = 0;}
       
  }
  return gg_flag;
}






void 
CLASS_NAME::begin_run( BelleEvent*, int *status ) {
  *status = 0;

  // read DB constants
  //  int nword = BsCouTab ( BELLE_RUNHEAD );
  if ( BsCouTab ( BELLE_RUNHEAD ) ) {
    Belle_runhead_Manager& runhead_mgr = Belle_runhead_Manager::get_manager();
    Belle_runhead& runhead(runhead_mgr[0]);
    ExpMC = runhead.ExpMC();
    ExpNo = runhead.ExpNo();
    RunNo = runhead.RunNo();
  }
  //
  BeamEnergy::begin_run();
  IpProfile::begin_run(); // Get IP profile data from $BELLE_POSTGRES_SERVER 
  IpProfile::dump();      // Dump IP profile data to STDOUT (optional) 
  eid::init_data();

  Hamlet::begin_run( (Hamlet::tag_method)18 );
  //
  // Get Belle moninal beam  energy 
  //
  int nword_bnb = BsCouTab ( BELLE_NOMINAL_BEAM );
  HepLorentzVector Beam_CM_BNB;
  Hep3Vector CMboost_BNB;
  if ( nword_bnb ) {
    Belle_nominal_beam_Manager& bnb_mgr = Belle_nominal_beam_Manager::get_manager();
    Belle_nominal_beam& bnb(bnb_mgr[0]);

    m_eler_bnb = -bnb.pz_low();
    m_eher_bnb =  bnb.pz_high();
    Beam_CM_BNB = HepLorentzVector (-m_eler_bnb*sin(theta),
				    0.,
				    -m_eher_bnb*cos(theta)+m_eler_bnb,
				    m_eher_bnb+m_eler_bnb);
    m_ebeam_bnb = Beam_CM_BNB.mag()/2.;
    //    CMboost_BNB = Beam_CM_BNB.boostVector();
  }
}






void 
CLASS_NAME::hist_def ( void ) {

  extern BelleTupleManager* BASF_Histogram;
  BelleTupleManager* tm = BASF_Histogram;
  //

  char *ntlist01;
  char *ntlist02;
  char *ntlist03;
  char *ntlist04;
  ntlist01 =
    "expno runno evtno "
    "t01 t02 t0n "
    "ebeam "
    "de mbc "

    "g1_lmom g1_cmom g1_lcos g1_ccos g1_e925 "
    "g2_lmom g2_cmom g2_lcos g2_ccos g2_e925 "
    "gg_lphi gg_cphi gg_lang gg_cang "

    "g1_t g2_t "

    "gam01_e gam02_e "
 
    "mc_b2gg mc_gg mc_g1 mc_g2 mc_g1m mc_g2m mc_g1gm mc_g2gm "
 
    "ncand "

    "gg_mass "

    "veto veto1 veto2 "
    //
    //
    "g1_pi0 g2_pi0 g1_eta g2_eta "
    "g4_pi0 g5_pi0 g4_eta g5_eta "
    "m_pi1 m_pi2 m_et1 m_et2 m_pi0 m_eta "
    "m_pi3 m_pi4 m_pi5 m_et3 m_et4 m_et5 "
    "prob_pi1 prob_pi2 prob_et1 prob_et2 prob_pi0 prob_eta "
    "clos_pi1 clos_pi2 clos_et1 clos_et2 clos_pi0 clos_eta "
    "clos_pi4 clos_pi5 clos_et4 clos_et5 clos_pi3 clos_et3 "

    "cosb dz qr "
    "k0mm2 k0et "
    "k0hoo0 k0hoo1 k0hoo2 k0hoo3 k0hoo4 "
    "k0hso00 k0hso01 k0hso02 k0hso03 k0hso04 "
    "k0hso10 k0hso12 k0hso14 k0hso20 k0hso22 k0hso24 "
    "k1mm2 k1et "
    "k1hoo0 k1hoo1 k1hoo2 k1hoo3 k1hoo4 "
    "k1hso00 k1hso01 k1hso02 k1hso03 k1hso04 "
    "k1hso10 k1hso12 k1hso14 k1hso20 k1hso22 k1hso24 ";

    //"lundgmo1 lundgmo2 ;

  ntlist02 =
    "expno runno evtno t01 t02 t0n ";

  ntlist03 =

    "veto veto1 veto2 "
    //
    //
    "g1_pi0 g2_pi0 g1_eta g2_eta "
    "g4_pi0 g5_pi0 g4_eta g5_eta "
    "m_pi1 m_pi2 m_et1 m_et2 m_pi0 m_eta "
    "m_pi3 m_pi4 m_pi5 m_et3 m_et4 m_et5 "
    "prob_pi1 prob_pi2 prob_et1 prob_et2 prob_pi0 prob_eta "
    "clos_pi1 clos_pi2 clos_et1 clos_et2 clos_pi0 clos_eta "
    "clos_pi4 clos_pi5 clos_et4 clos_et5 clos_pi3 clos_et3 "
    "mc_gg gg_mass gg_le gg_ce gg_lmom gg_cmom gg_lcos gg_ccos ";
    //"lundgmo1 lundgmo2" ;

  ntlist04 =
    "expno runno evtno "

    "de mbc "

    "veto veto1 veto2 "; 
    //gg_mass g1_e g2_e mc_gg gg_mass gg_le gg_ce gg_lmom gg_cmom gg_lcos gg_ccos ";

  T[0] = tm->ntuple( "b2gg",     ntlist01, 1);
  T[1] = tm->ntuple( "b2gg evt", ntlist02, 2);
  T[2] = tm->ntuple( "b2gg gg",  ntlist03, 3);
  T[3] = tm->ntuple( "b2gg gg2",  ntlist04, 4);

  H[0] =tm->histogram("Statistics",  10, 0.0, 10.0, 10);
  H[10] =tm->histogram("gg_mass",  100, 0.0, 10.0);
  H[11] =tm->histogram("zzz",  100, 0.0, 10);

}






int 
CLASS_NAME::rec_b2gg( std::vector<Particle> &list_b ){

  for( std::vector<Particle>::iterator iii = list_b.begin();
       iii != list_b.end(); ++iii ){

    HepLorentzVector p4lab_b = iii->p();
    HepLorentzVector p4cms_b = p4lab_b;
    p4cms_b.boost(CMBoost);

    double mbc2 = pow( ebeam, 2. ) - p4cms_b.vect().mag2();
    if (mbc2<0){
      list_b.erase(iii); 
      iii--;
      continue;
    }
    double mbc = sqrt( mbc2 );
    double de  = p4cms_b.t() - ebeam;
    if (fabs(de)>0.5 || mbc<5.20 ){ 
      list_b.erase(iii); 
      iii--;
    }

  }
  return list_b.size();
}





int
CLASS_NAME::u_mc_type( const Particle & ptcl_gam01,
		       const Particle & ptcl_gam02){


   Mdst_gamma gamma01;
   Mdst_gamma gamma02;

   int mc_b2gg = 0;

   gamma01 = ptcl_gam01.mdstGamma();
   gamma02 = ptcl_gam02.mdstGamma();

   const Mdst_gamma *ggam01  = &gamma01;
   const Mdst_gamma *ggam02  = &gamma02;

   mc_genid_cand[0] = mc_getid_gam(ptcl_gam01);
   mc_genid_cand[1] = mc_getid_gam(ptcl_gam02);
   mc_genid_cand[2] = mc_getid_mother_gam(ggam01);
   mc_genid_cand[3] = mc_getid_mother_gam(ggam02);
   mc_genid_cand[4] = mc_getid_gmother_gam(ggam01); 
   mc_genid_cand[5] = mc_getid_gmother_gam(ggam02);
   mc_genid_cand[6] = mc_getid_ggmother_gam(ggam01);
   mc_genid_cand[7] = mc_getid_ggmother_gam(ggam02);
   mc_genid_cand[8] = mc_getid_gggmother_gam(ggam01);
   mc_genid_cand[9] = mc_getid_gggmother_gam(ggam02);    
    


    if ( (mc_genid_cand[0]==mc_genid[0]) && (mc_genid_cand[1]==mc_genid[1]) ||
	  (mc_genid_cand[1]==mc_genid[0]) && (mc_genid_cand[0]==mc_genid[1]) ) {
	mc_b2gg = 1;
    }else if ( (mc_genid_cand[2]==mc_genid[0]) && (mc_genid_cand[1]==mc_genid[1]) ||
		(mc_genid_cand[1]==mc_genid[0]) && (mc_genid_cand[2]==mc_genid[1]) ||
		(mc_genid_cand[0]==mc_genid[0]) && (mc_genid_cand[3]==mc_genid[1]) ||
		(mc_genid_cand[3]==mc_genid[0]) && (mc_genid_cand[0]==mc_genid[1]) ||
		(mc_genid_cand[2]==mc_genid[0]) && (mc_genid_cand[3]==mc_genid[1]) ||
		(mc_genid_cand[3]==mc_genid[0]) && (mc_genid_cand[2]==mc_genid[1]) ) {
	mc_b2gg = 2;
    }else if ( (mc_genid_cand[4]==mc_genid[0]) && (mc_genid_cand[1]==mc_genid[1]) ||
		(mc_genid_cand[1]==mc_genid[0]) && (mc_genid_cand[4]==mc_genid[1]) ||
		(mc_genid_cand[0]==mc_genid[0]) && (mc_genid_cand[5]==mc_genid[1]) ||
		(mc_genid_cand[5]==mc_genid[0]) && (mc_genid_cand[0]==mc_genid[1]) ||
		(mc_genid_cand[4]==mc_genid[0]) && (mc_genid_cand[3]==mc_genid[1]) ||
		(mc_genid_cand[3]==mc_genid[0]) && (mc_genid_cand[4]==mc_genid[1]) ||
		(mc_genid_cand[2]==mc_genid[0]) && (mc_genid_cand[5]==mc_genid[1]) ||
		(mc_genid_cand[5]==mc_genid[0]) && (mc_genid_cand[2]==mc_genid[1]) ||
		(mc_genid_cand[4]==mc_genid[0]) && (mc_genid_cand[5]==mc_genid[1]) ||
		(mc_genid_cand[5]==mc_genid[0]) && (mc_genid_cand[4]==mc_genid[1]) ) {
	mc_b2gg = 3; 
    }else if ( (mc_genid_cand[6]==mc_genid[0]) && (mc_genid_cand[1]==mc_genid[1]) ||
		(mc_genid_cand[1]==mc_genid[0]) && (mc_genid_cand[6]==mc_genid[1]) ||
		(mc_genid_cand[0]==mc_genid[0]) && (mc_genid_cand[7]==mc_genid[1]) ||
		(mc_genid_cand[7]==mc_genid[0]) && (mc_genid_cand[0]==mc_genid[1]) ||
		(mc_genid_cand[6]==mc_genid[0]) && (mc_genid_cand[3]==mc_genid[1]) ||
		(mc_genid_cand[3]==mc_genid[0]) && (mc_genid_cand[6]==mc_genid[1]) ||
		(mc_genid_cand[2]==mc_genid[0]) && (mc_genid_cand[7]==mc_genid[1]) ||
		(mc_genid_cand[7]==mc_genid[0]) && (mc_genid_cand[2]==mc_genid[1]) ||
		(mc_genid_cand[6]==mc_genid[0]) && (mc_genid_cand[5]==mc_genid[1]) ||
		(mc_genid_cand[5]==mc_genid[0]) && (mc_genid_cand[6]==mc_genid[1]) ||
		(mc_genid_cand[4]==mc_genid[0]) && (mc_genid_cand[7]==mc_genid[1]) ||
		(mc_genid_cand[7]==mc_genid[0]) && (mc_genid_cand[4]==mc_genid[1]) ||
		(mc_genid_cand[6]==mc_genid[0]) && (mc_genid_cand[7]==mc_genid[1]) ||
		(mc_genid_cand[7]==mc_genid[0]) && (mc_genid_cand[6]==mc_genid[1]) ) {
	mc_b2gg = 4;
      }

/*
    std::vector<int> mc_genid_cand;
   
    mc_genid_cand.push_back(id_gam1);         //0
    mc_genid_cand.push_back(id_gam2);         //1
    mc_genid_cand.push_back(id_mother1);      //2
    mc_genid_cand.push_back(id_mother2);      //3
    mc_genid_cand.push_back(id_gmother_gam1); //4
    mc_genid_cand.push_back(id_gmother_gam2); //5

    if ( ((mc_genid_cand.at(0)==mc_genid[0]) && (mc_genid_cand.at(1)==mc_genid[1]) ||
	  (mc_genid_cand.at(0)==mc_genid[1]) && (mc_genid_cand.at(1)==mc_genid[0])) ) {
	mc_b2gg = 1;
    }else if ( ((mc_genid_cand.at(0)==mc_genid[0]) && (mc_genid_cand.at(3)==mc_genid[1]) ||
		(mc_genid_cand.at(3)==mc_genid[1]) && (mc_genid_cand.at(0)==mc_genid[0])) ) {
	mc_b2gg = 2;
    }else if ( ((mc_genid_cand.at(0)==mc_genid[0]) && (mc_genid_cand.at(5)==mc_genid[1]) ||
		(mc_genid_cand.at(5)==mc_genid[1]) && (mc_genid_cand.at(0)==mc_genid[0])) ) {
	mc_b2gg = 3;
    }else if ( ((mc_genid_cand.at(2)==mc_genid[0]) && (mc_genid_cand.at(1)==mc_genid[1]) ||
		(mc_genid_cand.at(1)==mc_genid[1]) && (mc_genid_cand.at(2)==mc_genid[0])) ) {
	mc_b2gg = 4;
    }else if ( ((mc_genid_cand.at(2)==mc_genid[0]) && (mc_genid_cand.at(3)==mc_genid[1]) ||
		(mc_genid_cand.at(3)==mc_genid[1]) && (mc_genid_cand.at(2)==mc_genid[0])) ) {
	mc_b2gg = 5;
    }else if ( ((mc_genid_cand.at(2)==mc_genid[0]) && (mc_genid_cand.at(5)==mc_genid[1]) ||
		(mc_genid_cand.at(5)==mc_genid[1]) && (mc_genid_cand.at(2)==mc_genid[0])) ) {
	mc_b2gg = 6;
    }else if ( ((mc_genid_cand.at(4)==mc_genid[0]) && (mc_genid_cand.at(1)==mc_genid[1]) ||
		(mc_genid_cand.at(1)==mc_genid[1]) && (mc_genid_cand.at(4)==mc_genid[0])) ) {
	mc_b2gg = 7;
    }else if ( ((mc_genid_cand.at(4)==mc_genid[0]) && (mc_genid_cand.at(3)==mc_genid[1]) ||
		(mc_genid_cand.at(3)==mc_genid[1]) && (mc_genid_cand.at(4)==mc_genid[0])) ) {
	mc_b2gg = 8;
    }else if ( ((mc_genid_cand.at(4)==mc_genid[0]) && (mc_genid_cand.at(5)==mc_genid[1]) ||
		(mc_genid_cand.at(5)==mc_genid[1]) && (mc_genid_cand.at(4)==mc_genid[0])) ) {
	mc_b2gg = 9;
    }else if ( (mc_genid_cand.at(0)==mc_genid[0]) || (mc_genid_cand.at(0)==mc_genid[1]) ) {
	mc_b2gg = 10;
    }else if ( (mc_genid_cand.at(1)==mc_genid[0]) || (mc_genid_cand.at(1)==mc_genid[1]) ) {
	mc_b2gg = 11;
      }
*/

  return mc_b2gg;
}







int
CLASS_NAME::mc_geninfo( void ){

  //
  // generated data 
  //
  int n_sig=0;

  Gen_hepevt_Manager &HepMgr = Gen_hepevt_Manager::get_manager();

  for (Gen_hepevt_Manager::iterator hi = HepMgr.begin();
       hi!= HepMgr.end(); hi++) {

    Gen_hepevt g_b = *hi;
    if ( abs(g_b.idhep()) != 511 ) { continue; }

    //
    // B0 -> gamma  gamma 
    //
    int n_dau = g_b.da(1)-g_b.da(0)+1;

    if ( n_dau != 2 ) { continue; }

    Gen_hepevt g_d01 = HepMgr(Panther_ID(g_b.da(0)));
    Gen_hepevt g_d02 = HepMgr(Panther_ID(g_b.da(1)));
    Gen_hepevt g_gam01;
    Gen_hepevt g_gam02;
    if ( g_d01.idhep()!=22 || g_d02.idhep()!=22 ){continue;}
    if (g_d01.E()>g_d02.E()){
      g_gam01 = g_d01;
      g_gam02 = g_d02;
    }else{
      g_gam01 = g_d02;
      g_gam02 = g_d01;
    }
    // fill ID 
    mc_genid[0] = g_gam01.get_ID(); 
    mc_genid[1] = g_gam02.get_ID(); 

    n_sig++;
  }
  return n_sig;
}
#if defined(BELLE_NAMESPACE)
} // namespace Belle 
#endif
